- Загрузите датасет life_expectancy_data.RDS (лежит в папке домашнего
задания). Это данные с основными показателями, через которые
высчитывается ожидаемая продолжительности жизни по метрике World
Development Indicator на уровне стран. В данных оставлены строки,
относящиеся к положению женщин в 2019 г.
data <- read_rds("C:/Users/79214/Documents/BioStat_2023/GitHub/life_expectancy_data.RDS")
summary(data)
## Country Year Gender Life expectancy
## Length:195 Min. :2019 Length:195 Min. :55.49
## Class :character 1st Qu.:2019 Class :character 1st Qu.:70.02
## Mode :character Median :2019 Mode :character Median :77.55
## Mean :2019 Mean :75.52
## 3rd Qu.:2019 3rd Qu.:80.95
## Max. :2019 Max. :88.10
## Unemployment Infant Mortality GDP GNI
## Min. : 0.178 Min. : 1.40 Min. :1.884e+08 Min. :3.754e+08
## 1st Qu.: 3.735 1st Qu.: 5.35 1st Qu.:1.117e+10 1st Qu.:1.094e+10
## Median : 5.960 Median :13.50 Median :3.967e+10 Median :4.009e+10
## Mean : 8.597 Mean :19.61 Mean :4.660e+11 Mean :4.864e+11
## 3rd Qu.:10.958 3rd Qu.:30.23 3rd Qu.:2.476e+11 3rd Qu.:2.457e+11
## Max. :36.442 Max. :75.80 Max. :2.143e+13 Max. :2.171e+13
## Clean fuels and cooking technologies Per Capita
## Min. : 0.00 Min. : 228.2
## 1st Qu.: 34.50 1st Qu.: 2165.3
## Median : 80.70 Median : 6624.8
## Mean : 65.98 Mean : 16821.0
## 3rd Qu.:100.00 3rd Qu.: 19439.7
## Max. :100.00 Max. :175813.9
## Mortality caused by road traffic injury Tuberculosis Incidence
## Min. : 0.00 Min. : 0.0
## 1st Qu.: 8.20 1st Qu.: 12.0
## Median :16.00 Median : 46.0
## Mean :17.06 Mean :103.8
## 3rd Qu.:24.00 3rd Qu.:138.5
## Max. :64.60 Max. :654.0
## DPT Immunization HepB3 Immunization Measles Immunization Hospital beds
## Min. :35.00 Min. :35.00 Min. :37.00 Min. : 0.200
## 1st Qu.:85.69 1st Qu.:81.31 1st Qu.:84.85 1st Qu.: 1.301
## Median :92.00 Median :91.00 Median :92.00 Median : 2.570
## Mean :87.99 Mean :86.76 Mean :87.31 Mean : 2.997
## 3rd Qu.:97.00 3rd Qu.:96.00 3rd Qu.:96.50 3rd Qu.: 3.773
## Max. :99.00 Max. :99.00 Max. :99.00 Max. :13.710
## Basic sanitation services Tuberculosis treatment Urban population
## Min. : 8.632 Min. : 0.00 Min. : 13.25
## 1st Qu.: 62.919 1st Qu.: 73.00 1st Qu.: 41.92
## Median : 91.144 Median : 82.00 Median : 58.76
## Mean : 77.380 Mean : 77.57 Mean : 59.12
## 3rd Qu.: 98.582 3rd Qu.: 88.00 3rd Qu.: 78.02
## Max. :100.000 Max. :100.00 Max. :100.00
## Rural population Non-communicable Mortality Sucide Rate continent
## Min. : 0.00 Min. : 4.40 Min. : 0.300 Africa :52
## 1st Qu.:21.98 1st Qu.:11.85 1st Qu.: 2.050 Americas:38
## Median :41.24 Median :17.20 Median : 3.500 Asia :42
## Mean :40.88 Mean :17.05 Mean : 4.802 Europe :48
## 3rd Qu.:58.08 3rd Qu.:22.10 3rd Qu.: 6.600 Oceania :15
## Max. :86.75 Max. :43.70 Max. :30.100
- Сделайте интерактивный plotly график любых двух нумерических
колонок. Раскрасть по колонке континента, на котором расположена
страна
plot_ly(
data = data[(data$Unemployment != 0) & (data$`Tuberculosis Incidence` != 0)],
x = ~ Unemployment,
y = ~ `Tuberculosis Incidence`,
color = ~ continent) %>%
layout(
title = 'Отношение уровня заболеваемости туберкулезом к уровню безработицы',
yaxis = list(title = 'Количество заболеваний туберкулезом',
zeroline = FALSE), # Уберём выделения нулевых осей по y
xaxis = list(title = 'Уровень безработицы',
zeroline = FALSE)) # Уберём выделения нулевых осей по y
- Проведите тест, на сравнение распределений колонки
Life expectancy между группами стран Африки и Америки. Вид
статистического теста определите самостоятельно. Визуализируйте
результат через библиотеку rstatix.
filter_data <- data %>%
filter(.$continent == 'Africa' | .$continent == 'Americas') %>%
select(Country, `Life expectancy`, continent)
filter_data %>%
ggqqplot(x = 'Life expectancy', facet.by = "continent")

#сравним распределение ожидаемой продолжительности жизни в Америке и Африке, для этого применим критерий Манна-Уитни
filter_data_aremicas <- filter_data %>%
filter(.$continent == 'Americas')
filter_data_africa <- filter_data %>%
filter(.$continent == 'Africa')
stat.test <- filter_data %>%
wilcox_test(`Life expectancy` ~ continent) %>%
add_xy_position(x = "continent")
stat.test #p-value << 0.05, значит можно отвергнуть нулевую гипотезу, следовательно ожидаемая продолжительность жизни разная в группах Африки и Америки
## # A tibble: 1 × 11
## .y. group1 group2 n1 n2 statistic p y.position groups xmin
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <name> <dbl>
## 1 Life exp… Africa Ameri… 52 38 107 6.34e-13 86.6 <chr> 1
## # ℹ 1 more variable: xmax <dbl>
ggboxplot(
filter_data,
x = "continent", y = 'Life expectancy' ,
ylab = 'Life expectancy', xlab = "Continent",
add = "jitter"
) +
labs(subtitle = get_test_label(stat.test, detailed = TRUE)) +
stat_pvalue_manual(stat.test, tip.length = 0)

- Сделайте новый датафрейм, в котором оставите все численные колонки
кроме
Year. Сделайте корреляционный анализ этих данных.
Постройте два любых типа графиков для визуализации корреляций.
new_data <- data %>%
select(-c(Country, Year, Gender, continent))
cor_new_data <- cor(new_data)
corrplot(cor_new_data, method = "number", type = "lower", tl.cex = 0.5, number.cex = 0.5)

ggpairs(new_data,
title = 'Correlations in Life expectancy dataset',progress = F) +
theme_minimal() +
scale_fill_manual(values = c('#69b3a2')) +
scale_colour_manual(values = c('#69b3a2'))

- Постройте иерархическую кластеризацию на этом датафрейме.
scaled_new_data <- scale(new_data)
scaled_clear_dist <- dist(scaled_new_data,
method = "euclidean")
as.matrix(scaled_clear_dist)[1:6,1:6]
## 1 2 3 4 5 6
## 1 0.000000 7.605708 6.331840 4.414874 6.645623 7.923487
## 2 7.605708 0.000000 2.624659 7.921597 3.357361 3.631018
## 3 6.331840 2.624659 0.000000 6.321666 4.350331 3.464837
## 4 4.414874 7.921597 6.321666 0.000000 8.095849 7.161240
## 5 6.645623 3.357361 4.350331 8.095849 0.000000 4.966244
## 6 7.923487 3.631018 3.464837 7.161240 4.966244 0.000000
data_clear_hc <- hclust(d = scaled_clear_dist,
method = "ward.D2")
fviz_dend(data_clear_hc,
cex = 0.1)
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

- Сделайте одновременный график heatmap и иерархической кластеризации.
Содержательно интерпретируйте результат.
pheatmap(scaled_new_data,
show_rownames = FALSE,
clustering_distance_rows = scaled_clear_dist,
clustering_method = "ward.D2",
cutree_rows = 5,
cutree_cols = length(colnames(scaled_new_data)),
angle_col = 45,
main = "Dendrograms for clustering rows and columns with heatmap")

#кластеризация помогает разбить колонки примерно на 5 групп, сами описательные признаки образуют 4 группы (т.е 4 региона), причем сильная корреляция по столбцам GDP и GNI в одном регионе. Между собой связны столбцы о безработице, количестве заболеваний туберкулезом, детская смертность, в друггую группы выделены столбцы 'Immunization', отдельная группа - количество мест в больницах и количество суицидов. По первой группе есть корреляция по одному региону.
- Проведите PCA анализ на этих данных. Проинтерпретируйте
результат.
data_full.pca <- prcomp(new_data,
scale = T)
summary(data_full.pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion 0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion 0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
## PC15 PC16 PC17 PC18 PC19
## Standard deviation 0.34546 0.26941 0.20224 0.06968 1.017e-15
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.000e+00
## Cumulative Proportion 0.99377 0.99759 0.99974 1.00000 1.000e+00
fviz_eig(data_full.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(data_full.pca, col.var = "contrib")

fviz_pca_var(data_full.pca,
select.var = list(contrib = 3),
col.var = "contrib") #очень сильная корреляция переменной 'Life expectations' с 1 компонентой

fviz_contrib(data_full.pca, choice = "var", axes = 1, top = 24)

fviz_contrib(data_full.pca, choice = "var", axes = 2, top = 24)

fviz_contrib(data_full.pca, choice = "var", axes = 3, top = 24)

#первые две компоненты объясняют 50% дисперсии компонент
#большой вклад в анализируемые данные компоненты вносят переменные: Unempoyment, sucide rate, tuberculosis treatment, hospital beds и т.д
#выделяются группы: "Immunization", 1 четверть и граница 2 и 3 четверти.
- Постройте biplot график для PCA. Раскрасьте его по значениям
континентов. Переведите его в
plotly. Желательно, чтобы при
наведении на точку, вы могли видеть название страны.
ggbiplot(data_full.pca,
scale=0, alpha = 0.1) +
theme_minimal()

data_clear_with_ch <- data %>%
select(-c(Country, Year, Gender))
ggplotly(ggbiplot(data_full.pca,
scale=0,
groups = as.factor(data_clear_with_ch$continent),
ellipse = T,
alpha = 0.2) +
theme_minimal())
- Дайте содержательную интерпретацию PCA анализу.
#дополнения к предыдущим объяснениям. Данные могу быть объяснены 3 переменными. Есть переменные, которые имеют отрицательную корреляцию. На последнем графике видны выбросы, данные можно кластеризовать по континентам, но это происходит не очень эффективно.
- Сравните результаты отображения точек между алгоритмами PCA и
UMAP.
umap_prep <- recipe(~., data = new_data) %>%
step_normalize(all_predictors()) %>%
step_umap(all_predictors()) %>%
prep() %>%
juice()
umap_prep %>%
ggplot(aes(UMAP1, UMAP2)) +
geom_point(aes(color = as.character(data_clear_with_ch$continent)),
alpha = 0.7, size = 2) +
labs(color = NULL)

#есть сходства в результатах отображения точек: точки Африканского континета аггрегируются в левом верхнем углу, как и в PCA, видна аггрегация Европейского континента (хотя оно не такое "сконцентрированное" как в PCA), нет таких сильных выбросов в UMAP. В UMAP точки находятся более плотно, лучше видны кластеры.
- Давайте самостоятельно увидим, что снижение размерности – это группа
методов, славящаяся своей неустойчивостью. Удалите 5 случайных колонок.
Проведите PCA анализ. Повторите результат 3 раза. Наблюдаете ли вы
изменения в куммулятивном проценте объяснённой вариации? В итоговом
представлении данных на биплотах? С чем связаны изменения между тремя
PCA?
#1
rand_1 <- sample(1:19, 5)
new_data_1 <- subset(new_data, select = -c(rand_1))
data_full.pca_1 <- prcomp(new_data_1,
scale = T)
summary(data_full.pca_1)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.3501 1.3777 1.15047 1.05237 0.97405 0.88640 0.79875
## Proportion of Variance 0.3945 0.1356 0.09454 0.07911 0.06777 0.05612 0.04557
## Cumulative Proportion 0.3945 0.5301 0.62460 0.70370 0.77147 0.82759 0.87317
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.6919 0.63884 0.55333 0.54901 0.3608 0.33158 0.20260
## Proportion of Variance 0.0342 0.02915 0.02187 0.02153 0.0093 0.00785 0.00293
## Cumulative Proportion 0.9074 0.93652 0.95838 0.97991 0.9892 0.99707 1.00000
fviz_eig(data_full.pca_1, addlabels = T, ylim = c(0, 40))

fviz_pca_var(data_full.pca_1, col.var = "contrib")

fviz_pca_var(data_full.pca_1,
select.var = list(contrib = 3),
col.var = "contrib") #объяснено 55,5% объясненной вариации

graph_1 <- ggplotly(ggbiplot(data_full.pca_1,
scale=0,
groups = as.factor(data_clear_with_ch$continent),
ellipse = T,
alpha = 0.2) +
theme_minimal())
#2
rand_2 <- sample(1:19, 5)
new_data_2 <- subset(new_data, select = -c(rand_2))
data_full.pca_2 <- prcomp(new_data_2,
scale = T)
summary(data_full.pca_2)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.4795 1.2080 1.14397 1.01342 0.96305 0.90564 0.75882
## Proportion of Variance 0.4391 0.1042 0.09348 0.07336 0.06625 0.05858 0.04113
## Cumulative Proportion 0.4391 0.5434 0.63683 0.71019 0.77644 0.83502 0.87615
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.68138 0.61773 0.58059 0.47287 0.3871 0.32157 0.27205
## Proportion of Variance 0.03316 0.02726 0.02408 0.01597 0.0107 0.00739 0.00529
## Cumulative Proportion 0.90932 0.93657 0.96065 0.97662 0.9873 0.99471 1.00000
fviz_eig(data_full.pca_2, addlabels = T, ylim = c(0, 40))

fviz_pca_var(data_full.pca_2, col.var = "contrib")

fviz_pca_var(data_full.pca_2,
select.var = list(contrib = 3),
col.var = "contrib") #объяснено 54,1% объясненной вариации

graph_2 <- ggplotly(ggbiplot(data_full.pca_2,
scale=0,
groups = as.factor(data_clear_with_ch$continent),
ellipse = T,
alpha = 0.2) +
theme_minimal())
#3
rand_3 <- sample(1:19, 5)
new_data_3 <- subset(new_data, select = -c(rand_3))
data_full.pca_3 <- prcomp(new_data_3,
scale = T)
summary(data_full.pca_3)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.3060 1.4008 1.3197 1.12170 1.05576 0.86576 0.76065
## Proportion of Variance 0.3798 0.1402 0.1244 0.08987 0.07962 0.05354 0.04133
## Cumulative Proportion 0.3798 0.5200 0.6444 0.73427 0.81389 0.86743 0.90876
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.69631 0.63187 0.50350 0.3062 0.20290 0.07002 3.191e-16
## Proportion of Variance 0.03463 0.02852 0.01811 0.0067 0.00294 0.00035 0.000e+00
## Cumulative Proportion 0.94339 0.97191 0.99001 0.9967 0.99965 1.00000 1.000e+00
fviz_eig(data_full.pca_3, addlabels = T, ylim = c(0, 40))

fviz_pca_var(data_full.pca_3, col.var = "contrib")

fviz_pca_var(data_full.pca_3,
select.var = list(contrib = 3),
col.var = "contrib") #объяснено 54,6% объясненной вариации

graph_3 <- ggplotly(ggbiplot(data_full.pca_3,
scale=0,
groups = as.factor(data_clear_with_ch$continent),
ellipse = T,
alpha = 0.2) +
theme_minimal())
graph_1
graph_2
graph_3
#при удалении любых 5 столбцов кумулятивный процент объясненной вариации примерно равно 55%, следовательно каждый удаленный столбец поднимает объясненную вариацию примерно на 1%. Эллипсы на биплотах стали более вытяннутыми, т.е данные стали более разбросаны. В зависимости от того, какие столбцы были удалены, расположение эллипсов отличается: эллипсы превращаются в круги, происходят перемещения эллипсов по графику. Все 3 графика достаточно сильно отличаются друг от друга.
- Давайте самостоятельно увидим, что снижение размерности – это группа
методов, славящаяся своей неустойчивостью. Создайте две дамми-колонки о
том: (1) принадлежит ли страна к африканскому континенту, (2) Океании.
Проведите PCA вместе с ними, постройте биплоты. Проинтерпрейтируйте
результат. Объясните, почему добавление дамми-колонок не совсем
корректно в случае PCA нашего типа.
dammi_data <- data %>%
mutate(Is_African = case_when(.$continent == 'Africa' ~ 1,
.$continent != 'Africa' ~ 0),
Is_Oceanian = case_when(.$continent == 'Oceania' ~ 1,
.$continent != 'Oceania' ~ 0))
dammi_data <- dammi_data %>%
select(-c(Country, Year, Gender, continent))
dammi_data.pca <- prcomp(dammi_data,
scale = T)
summary(dammi_data.pca)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 2.8400 1.4851 1.3960 1.28878 1.11384 1.01306 0.9514
## Proportion of Variance 0.3841 0.1050 0.0928 0.07909 0.05908 0.04887 0.0431
## Cumulative Proportion 0.3841 0.4891 0.5819 0.66098 0.72006 0.76893 0.8120
## PC8 PC9 PC10 PC11 PC12 PC13 PC14
## Standard deviation 0.92891 0.84328 0.70113 0.67720 0.62108 0.55110 0.46766
## Proportion of Variance 0.04109 0.03386 0.02341 0.02184 0.01837 0.01446 0.01041
## Cumulative Proportion 0.85312 0.88699 0.91040 0.93223 0.95060 0.96506 0.97548
## PC15 PC16 PC17 PC18 PC19 PC20
## Standard deviation 0.39263 0.35858 0.34148 0.26517 0.20133 0.06900
## Proportion of Variance 0.00734 0.00612 0.00555 0.00335 0.00193 0.00023
## Cumulative Proportion 0.98282 0.98894 0.99449 0.99784 0.99977 1.00000
## PC21
## Standard deviation 9.958e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion 1.000e+00
fviz_eig(dammi_data.pca, addlabels = T, ylim = c(0, 40))

fviz_pca_var(dammi_data.pca, col.var = "contrib")

fviz_pca_var(dammi_data.pca,
select.var = list(contrib = 3),
col.var = "contrib") #объяснено 55,5% объясненной вариации

ggplotly(ggbiplot(dammi_data.pca,
scale=0,
groups = as.factor(data_clear_with_ch$continent),
ellipse = T,
alpha = 0.2) +
theme_minimal())
#добавление двух дамми-колонок уменьшает кумулятивный процент объясненной вариации примерно на 2%, но точки остались на прежних местах, эллипсы остались такими же. То есть наши данные просто стали хуже объясняться глаными компонентами. Это происходит из-за того, что добавленные столбцы бинарные и не несут количественного и качественного смысла без категориальных столцов.